# test the sampling method
# 
# Author: guochun
###############################################################################

source("./sourceAll.R")
#source("./objects/sourceAll.R")
#source("./sampling/sourceAll.R")
library(spatstat)
library(FNN)
library(multicore)

data=getData("bci5.txt")
availible=data$status=="A" & data$sp == "ACALDI" & !is.na(data$gx) & !is.na(data$gy)
subdata=data[availible,]

#population=new("population",x=subdata$gx,y=subdata$gy,plotdim=c(1000,500))
population=new("population",x=sample(1:1000,100),y=sample(1:1000,100),plotdim=c(1000,1000))
#population=new("population",x=rep(seq(1,1000,50),20),y=rep(seq(1,1000,50),each=20),plotdim=c(1000,1000))
#the general one

source("./sourceAll.R")
angleorder=new("angleOrder",n=10,q=3,k=4,visual=T,time=100,overlap=FALSE)
sampling(angleorder,population,quick=T)

angleorder=new("angleOrder",n=10,q=3,k=4,visual=T,time=100,edge=80)
sampling(angleorder,population,quick=T)

angleorder=new("angleOrder",n=50,q=4,k=1,focusType="event",visual=T)
angleorder=new("angleOrder",n=50,q=4,k=1,focusType="point",visual=T)
sampling(angleorder,population,quick=T)
angleorder=new("angleOrder",n=72,q=1,k=1,visual=T,composite=T)
sampling(angleorder,population,quick=T)

angleorder=new("angleOrder",n=100,q=1,k=4)
sampling(angleorder,population,quick=T)
angleorder=new("angleOrder",n=10,q=2,k=1)
sampling(angleorder,population,quick=T)
angleorder=new("angleOrder",n=10,q=3,k=1)

angleorder=new("angleOrder",n=10,q=2,k=2)
sampling(angleorder,population,quick=T)

angleorder=new("angleOrder",n=100,q=4,k=3)
sampling(angleorder,population,quick=T)
re=list()
for (q in 1:4){
	for (k in 1:3){
		for(n in 100:200){
			#angleorder=new("angleOrder",n=n,q=q,k=k)
			angleorder=new("angleOrder",n=n,q=q,k=k,focusType="event")
			#for(i in 1:10){
				re[[k+(q-1)*3]]=sampling(angleorder,population,quick=T)
			#}
				
		}
	}
}


#can overlap
eventToEvent=new("eToe",time=100)
sampling(eventToEvent,population,T)
pointToEvent=new("pToe",time=100)
sampling(pointToEvent,population,T)
tSquareSample=new("tSquare",time=50)
sampling(tSquareSample,population,T)
pQuaterSample=new("pQuarter",time=10)
sampling(pQuaterSample,population,T)

#none overlap
eventToEvent=new("eToe",time=100,overlap=F)
sampling(eventToEvent,population,T)
pointToEvent=new("pToe",time=100,overlap=F)
sampling(pointToEvent,population,T)
tSquareSample=new("tSquare",time=50,overlap=F)
sampling(tSquareSample,population,T)
pQuaterSample=new("pQuarter",time=10,overlap=F)
sampling(pQuaterSample,population,T)

#border edge correction
eventToEvent=new("eToe",time=100,overlap=F,edge=30)
sampling(eventToEvent,population,T)
pointToEvent=new("pToe",time=100,overlap=F,edge=30)
sampling(pointToEvent,population,T)
tSquareSample=new("tSquare",time=50,overlap=F,edge=30)
sampling(tSquareSample,population,T)
pQuaterSample=new("pQuarter",time=10,overlap=F,edge=20)
sampling(pQuaterSample,population,T)

#test the fast calculation
sampling(eventToEvent,population,quick=T)
sampling(pointToEvent,population,quick=T)
tSquareSample=new("tSquare",time=100,overlap=F,edge=30)
sampling(tSquareSample,population,quick=T)
pQuaterSample=new("pQuarter",time=100,overlap=F,edge=20)
sampling(pQuaterSample,population,quick=T)
pQuaterSample2=new("pQuarter",time=100,overlap=F,edge=20,transectLine=FALSE)
sampling(pQuaterSample2,population,quick=T)
pointToEvent2=new("pToe",time=50,overlap=F,edge=30,transectLine=TRUE)
sampling(pointToEvent2,population,quick=T)
pointToEvent=new("pToe",time=50,overlap=F,edge=30,transectLine=TRUE,nonparametric=TRUE)

#test the kth nearest neighborhood distance
eventToEvent=new("eToe",time=100,overlap=F,edge=10,k=2)
sampling(eventToEvent,population,quick=T)
pointToEvent=new("pToe",time=100,overlap=F,edge=10,k=2)
sampling(pointToEvent,population,quick=T)
tSquareSample=new("tSquare",time=100,overlap=F,edge=10,k=2)
sampling(tSquareSample,population,quick=T)
pQuaterSample=new("pQuarter",time=30,overlap=F,edge=20,k=2)
sampling(pQuaterSample,population,quick=T)

#test the composite index
tSquareSample=new("tSquare",time=100,overlap=F,edge=10,k=1,composite=TRUE)
sampling(tSquareSample,population,quick=T)

source("./sourceAll.R")
#using random pattern to test the correction of our R code
eventToEvent=new("eToe",time=100,overlap=F,edge=10,k=1)
pointToEvent=new("pToe",time=100,overlap=F,edge=10,k=2)
tSquareSample=new("tSquare",time=100,overlap=F,edge=10,k=2)
pQuaterSample=new("pQuarter",time=25,overlap=F,edge=10,k=2)
diggleCpSample=new("diggleCp",time=100,overlap=F,edge=10,k=1)
tSquareSample=new("tSquare",time=100,overlap=F,edge=10,k=1,composite=TRUE)

#sampleMethod=eventToEvent
#sampleMethod=pointToEvent
#sampleMethod=tSquareSample
#sampleMethod=pQuaterSample
#sampleMethod=diggleCpSample
sampleMethod=new("angleOrder",n=10,q=2,k=2,time=100)
reptime=1000
#re=numeric()
#for(i in 1:reptime){
#np=rpois(1,1000)
ppo=ppp(x=runif(10000,0,1000),y=runif(10000,0,1000),window=owin(c(0,1000),c(0,1000)))
ppo=unique(ppo)
population=new("population",x=ppo$x,y=ppo$y,plotdim=c(1000,1000))
#population=new("population",x=runif(1000,0,1000),y=runif(1000,0,500),plotdim=c(1000,500))
re=unlist(lapply(1:reptime,function(x) pSizeTor(sampling(sampleMethod,population,quick=T))))
#re=c(re,sampling(sampleMethod,population,quick=T))
#}
mean(re)

re=numeric()
nsamples=seq(10,500,30)
for(j in 1:length(nsamples)){
	sampleMethod=new("diggleCp",time=nsamples[j],overlap=F,edge=10,k=1)
	re[j]=mean(unlist(lapply(1:reptime,function(x) pSizeTor(sampling(sampleMethod,population,quick=T)))))
}



